home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / IntermediatC / exponential sim #10 / exponential.c next >
Encoding:
C/C++ Source or Header  |  1990-11-15  |  6.8 KB  |  188 lines
  1. /*
  2.  *
  3.  * Nondeterministic Systems Computer Simulation 2
  4.  *
  5.  * Exponential Distribution
  6.  *
  7.  * This program simulates an exponential distribution.  It generates a number by the
  8.  * equation y = -ln(x), where x is a random number greater than 0 and less than or
  9.  * equal to 1 (since ln(0) is not defined).  The program also generates the theoretical
  10.  * frequency by obtaining the theoretical probability between two points and multiplying
  11.  * by the number of samples.  From the frequency, the distribution function can be
  12.  * obtained by summing the frequencies from 0 to the current point and dividing by the
  13.  * total number of samples taken.
  14.  *
  15.  * Develop a program that performs the function y = -ln(x) where x is a random
  16.  * real number between 0 and 1 (but it cannot equal 0, although 1 is allowed).
  17.  * Keep track of the number of times a value y falls between an arbitrary interval
  18.  * (which you decide upon).  This is the probability that the random variable y
  19.  * falls within the interval.  The distribution function is the probability that
  20.  * y <= a specified value.  This can be accomplished by summing all densities
  21.  * (what was calculated previously) from the beginning up to the current y.
  22.  *
  23.  * The random variable y should be calculated thousands of times and the table
  24.  * of distributions should be printed to a file.
  25.  *
  26.  */
  27.  
  28.  
  29. /* ------------------------------------------------------------------------------- */
  30.  
  31.  
  32. #include <stdio.h>
  33. #include <stdlib.h>
  34. #include <math.h>
  35. #include <time.h>
  36.  
  37.  
  38. /* ------------------------------------------------------------------------------- */
  39.  
  40.  
  41. #define MAX_ITERATIONS    500000L        /* number of samples generated */
  42. #define MAX_INTERVALS    100            /* maximum number of samples intervals */
  43. #define DISCRETE        10.0        /* number of discrete samples for an interval of 1:
  44.                                         10 samples between 0 and 1, another 10 between
  45.                                         1 and 2, etc. */
  46. #define ONE_TENTH        0            /* indexes for probability array */
  47. #define ONE_HALF        1
  48. #define ONE                2
  49. #define TWO                3
  50.  
  51.  
  52. /* ------------------------------------------------------------------------------- */
  53.  
  54.  
  55. main()
  56. {
  57.     long i;                                /* loop counter */
  58.     long count;                            /* variables to sum the frequencies */
  59.     long count2;
  60.     double y;                            /* random variable */
  61.     double divisor = RAND_MAX + 1.0;    /* used to divide random number; done as variable
  62.                                             to save addition and double conversion time */
  63.     long density[MAX_INTERVALS];        /* array stores sampled density */
  64.     long probability[TWO+1];            /* array to store probability frequencies */
  65.     long theoretical[MAX_INTERVALS];    /* array to store theoretical sampling */
  66.     double p_theoretical[TWO+1];        /* array of theoretical probabilities */
  67.     FILE *fp;                            /* file control block pointer */
  68.     
  69.     srand(time(NULL));                    /* seeds the random number generator */
  70.     
  71. /*
  72.     Clear arrays that will be incremented.
  73. */
  74.  
  75.     for (i=0L; i<(long) MAX_INTERVALS; i++)
  76.         density[i] = 0;
  77.     for (i=0L; i<(TWO + 1L); i++)
  78.         probability[i] = 0;
  79.  
  80. /*
  81.     For comparison, the theoretical probability is generated below.  This probability
  82.     is applied to an interval by multiplying it with the maximum iterations, thus a
  83.     number of theoretical samples are created to compare with the simulation.
  84. */
  85.  
  86.     for (i=0L; i<(long) MAX_INTERVALS; i++)
  87.     {
  88.         y = -( exp( -(i+1)/DISCRETE ) - exp( -i/DISCRETE ) );
  89.         theoretical[i] = (long) floor(y * MAX_ITERATIONS + 0.5);
  90.     }
  91.     p_theoretical[ONE_TENTH] = -( exp( -0.1 ) - exp( 0.0 ) );
  92.     p_theoretical[ONE_HALF] = -( exp( -0.5 ) - exp( 0.0 ) );
  93.     p_theoretical[ONE] = -( exp( -1.0 ) - exp( 0.0 ) );
  94.     p_theoretical[TWO] = -( exp( -2.0 ) - exp( 0.0 ) );
  95.     
  96. /*
  97.     The code below does the actual simulation.  A random number is generated by use
  98.     of rand() between 0 and RAND_MAX (2^31 - 1).  This number must be converted to
  99.     to a number varying between 0 and 1 (but cannot equal 0).  Thus, rand() is divided
  100.     by (RAND_MAX + 1) and subtracted from 1.  If the value is less than certain
  101.     numbers, their probability is incremented.
  102. */
  103.  
  104.     for (i=0L; i<MAX_ITERATIONS; i++)
  105.     {
  106.         y = -log(1.0 - ((double) rand() / divisor));
  107.         if (y <= 0.1)
  108.             probability[ONE_TENTH]++;
  109.         if (y <= 0.5)
  110.             probability[ONE_HALF]++;
  111.         if (y <= 1.0)
  112.             probability[ONE]++;
  113.         if (y <= 2.0)
  114.             probability[TWO]++;
  115.             
  116. /*
  117.     Now, the random variable is multiplied by DISCRETE.  This allows for sampling.  A
  118.     sample is taken DISCRETE times between two successive integers (such as 0 and 1 or 
  119.     1 and 2, etc.).  If the random variable is within the sample range (less than
  120.     MAX_INTERVALS), the density is incremented.  The sampling needs to be done
  121.     since the function provided is continuous, but a this computer simulation only deals
  122.     with discrete numbers.
  123. */
  124.  
  125.         y *= DISCRETE;
  126.         if ( ((int) y) < MAX_INTERVALS )
  127.             density[(int) y]++;
  128.     }
  129.     
  130. /*
  131.     The data is then output to a file.  Three output files are created.  Two will just
  132.     contain the raw data with no labels.  This data will be copied and placed into Excel.
  133.     The other data is labeled so that reading is quite easy.  This file is included
  134.     with the report.
  135. */
  136.       
  137.     fp = fopen("Sim density.dat", "w");                /* raw data for simulation density */
  138.     for (i=0; i<(long) MAX_INTERVALS; i++)
  139.         fprintf(fp, "%ld\n", density[i]);
  140.     fclose(fp);
  141.  
  142.     fp = fopen("Theor density.dat", "w");            /* raw data for theoretical density */
  143.     for (i=0; i<(long) MAX_INTERVALS; i++)
  144.         fprintf(fp, "%ld\n", theoretical[i]);
  145.     fclose(fp);
  146.     
  147.     count = 0L;
  148.     fp = fopen("Sim distribution.dat", "w");    /* raw data for simulation distribution */
  149.     for (i=0; i<(long) MAX_INTERVALS; i++)
  150.     {
  151.         count += density[i];
  152.         fprintf(fp, "%lf\n", (double) count / (double) MAX_ITERATIONS);
  153.     }
  154.     fclose(fp);
  155.  
  156.     count = 0L;
  157.     fp = fopen("Theor distribution.dat", "w");    ./* raw data for theoretical distribution */
  158.     for (i=0; i<(long) MAX_INTERVALS; i++)
  159.     {
  160.         count += theoretical[i];
  161.         fprintf(fp, "%lf\n", (double) count / (double) MAX_ITERATIONS);
  162.     }
  163.     fclose(fp);
  164.  
  165.     fp = fopen("Exponential.out", "w");
  166.     fprintf(fp, "Minimum Value   Maximum Value   Samples   Theor Samples");
  167.     fprintf(fp, "   Distribution   Theor Distribution\n");
  168.     count = 0L;
  169.     count2 = 0L;
  170.     for (i=0; i<(long) MAX_INTERVALS; i++)
  171.     {
  172.         count += density[i];
  173.         count2 += theoretical[i];
  174.         fprintf(fp, "%8.1lf %16.1lf %12ld %12ld %16.6lf %17.6lf\n", i/DISCRETE, 
  175.             (i+1)/DISCRETE, density[i], theoretical[i], (double) count / 
  176.             (double) MAX_ITERATIONS, (double) count2 / (double) MAX_ITERATIONS);
  177.     }
  178.     fprintf(fp, "\n\nP(y<=0.1) = %lf     P(theoretical)(y<=0.1) = %lf\n", 
  179.         (double) probability[ONE_TENTH]/(double) MAX_ITERATIONS, p_theoretical[ONE_TENTH]);     
  180.     fprintf(fp, "P(y<=0.5) = %lf     P(theoretical)(y<=0.5) = %lf\n", 
  181.         (double) probability[ONE_HALF]/(double) MAX_ITERATIONS, p_theoretical[ONE_HALF]);     
  182.     fprintf(fp, "P(y<=1.0) = %lf     P(theoretical)(y<=1.0) = %lf\n", 
  183.         (double) probability[ONE]/(double) MAX_ITERATIONS, p_theoretical[ONE]);     
  184.     fprintf(fp, "P(y<=2.0) = %lf     P(theoretical)(y<=2.0) = %lf\n", 
  185.         (double) probability[TWO]/(double) MAX_ITERATIONS, p_theoretical[TWO]);     
  186.     fclose(fp);
  187. }
  188.